Skrypty geoprzetwarzania

Klasteryzacja danych

Author

Krzysztof Dyba

library("terra")
library("rstac")
options(timeout = 600)

Pozyskanie danych

Niniejszą analizę wykonamy na podstawie danych z Sentinela-2. Jako źródło danych wykorzystamy katalog STAC i usługę Earth Search. Przeprowadzimy również proste filtrowanie dostępnych produktów definiując w zapytaniu: kolekcję (collections), zakres przestrzenny w układzie WGS 84 (bbox), interwał czasowy w standardzie RFC 3339 (datetime) oraz atrybut zachmurzenia sceny (eo:cloud_cover).

stac_source = stac("https://earth-search.aws.element84.com/v1")
stac_source |>
  stac_search(
    collections = "sentinel-2-c1-l2a",
    bbox = c(22.5, 51.1, 22.6, 51.2),
    datetime = "2023-03-01T00:00:00Z/2023-10-31T00:00:00Z") |>
  ext_query(`eo:cloud_cover` < 10) |>
  post_request() -> obrazy

Jako wynik zapytania zostały zwrócone metadane scen spełniające zadane warunki. Następnie możemy wybrać przykładową scenę z 21 września 2023 r. oraz ograniczyć liczbę kanałów do czterech podstawowych w rozdzielczości 10 m, tj. niebieskiego, zielonego, czerwonego i bliskiej podczerwieni w celu uproszczenia analizy. Finalnie, używając funkcji assets_url() zostaną zwrócone odpowiednie odnośniki do pobrania rastrów.

kanaly = c("blue", "green", "red", "nir")
obrazy |>
  items_filter(properties$`s2:tile_id` == "S2A_OPER_MSI_L2A_TL_2APS_20230921T151458_A043076_T34UEB_N05.09") |>
  assets_select(asset_names = kanaly) |>
  assets_url() -> sentinel
sentinel
[1] "https://e84-earth-search-sentinel-data.s3.us-west-2.amazonaws.com/sentinel-2-c1-l2a/34/U/EB/2023/9/S2A_T34UEB_20230921T094329_L2A/B02.tif"
[2] "https://e84-earth-search-sentinel-data.s3.us-west-2.amazonaws.com/sentinel-2-c1-l2a/34/U/EB/2023/9/S2A_T34UEB_20230921T094329_L2A/B03.tif"
[3] "https://e84-earth-search-sentinel-data.s3.us-west-2.amazonaws.com/sentinel-2-c1-l2a/34/U/EB/2023/9/S2A_T34UEB_20230921T094329_L2A/B08.tif"
[4] "https://e84-earth-search-sentinel-data.s3.us-west-2.amazonaws.com/sentinel-2-c1-l2a/34/U/EB/2023/9/S2A_T34UEB_20230921T094329_L2A/B04.tif"

Zauważ, że kanał 8 (bliska podczerwień) jest przed kanałem 4 (czerwony). Należy mieć na uwadze, że nieodpowiednia kolejność może spowodować problemy w dalszej części projektu.

W kolejnym kroku, używając funkcji dir.create() stworzymy nowy katalog na dysku, do którego zostaną pobrane zobrazowania. Oprócz tego, należy zdefiniować ścieżki i nazwy plików. W tym celu będą przydatne dwie funkcje – file.path() do stworzenia ścieżek do plików oraz basename() do wyodrębnienia nazwy pliku wraz z rozszerzeniem z URL.

dir.create("sentinel")
rastry = file.path("sentinel", basename(sentinel))

Teraz możemy pobrać nasze dane używając funkcji download.file() w pętli.

for (i in seq_along(sentinel)) {
  download.file(sentinel[i], rastry[i], mode = "wb")
}

Przygotowanie danych

Po pobraniu danych musimy stworzyć listę plików (rastrów), które zamierzamy wczytać. W tym celu możemy wykorzystać funkcję list.files(), która jako argument przyjmuje ścieżkę do folderu z plikami. Oprócz tego musimy wskazać jaki rodzaj plików chcemy wczytać (pattern = "\\.tif$") oraz zwrócić pełne ścieżki do plików (full.names = TRUE).

rastry = list.files("sentinel", pattern = "\\.tif$", full.names = TRUE)
rastry
[1] "sentinel/B02.tif" "sentinel/B03.tif" "sentinel/B04.tif" "sentinel/B08.tif"

Kiedy utworzyliśmy już listę plików, możemy je wczytać za pomocą funkcji rast(). Pobrane zobrazowania pokrywają duży obszar (około 12 000 km\(^2\)). Dla uproszczenia analizy zdefiniujmy mniejszy zakres przestrzenny do wczytania używając funkcji ext() oraz przekazując obiekt SpatExtent jako argument win w funkcji rast().

Zwróć uwagę, że do zdefiniowania zakresu przestrzennego podczas wyszukiwania danych użyliśmy układu WGS 84, natomiast teraz wymagany jest rzeczywisty układ rastra, tj. EPSG:32634. Można to sprawdzić w metadanych katalogu STAC (obiekt obrazy) lub używając funkcji describe() (wymaga ścieżki do pliku).

bbox = ext(510000, 540000, 5630000, 5650000)
r = rast(rastry, win = bbox)

Możemy również zmienić nazwy kanałów spektralnych. Przed tą operacją należy się upewnić czy kanały zostały wczytane w prawidłowej kolejności.

names(r) = kanaly
r
class       : SpatRaster 
dimensions  : 2000, 3000, 4  (nrow, ncol, nlyr)
resolution  : 10, 10  (x, y)
window      : 510000, 540000, 5630000, 5650000  (xmin, xmax, ymin, ymax)
coord. ref. : WGS 84 / UTM zone 34N (EPSG:32634) 
sources     : B02.tif  
              B03.tif  
              B04.tif  
              B08.tif  
names       :   blue,  green,    red,    nir 
min values  : >    1, >  352, >  677, >  910 
max values  : 18560<, 17648<, 17056<, 16514< 

Następnie możemy sprawdzić podstawowe statystyki opisowe używając funkcji summary(). W przypadku dużych rastrów, zostanie automatycznie wykorzystana próba 100 tys. komórek.

summary(r)
Warning: [summary] used a sample
      blue              green              red              nir        
 Min.   :-0.00080   Min.   :0.00440   Min.   :0.0027   Min.   :0.0008  
 1st Qu.: 0.01510   1st Qu.:0.03260   1st Qu.:0.0211   1st Qu.:0.1951  
 Median : 0.03130   Median :0.06070   Median :0.0497   Median :0.2371  
 Mean   : 0.03579   Mean   :0.06165   Mean   :0.0622   Mean   :0.2491  
 3rd Qu.: 0.05170   3rd Qu.:0.08420   3rd Qu.:0.0958   3rd Qu.:0.2928  
 Max.   : 0.49840   Max.   :0.49280   Max.   :0.5080   Max.   :0.6909  

Zasadniczo, wartości odbicia spektralnego mieszczą się w przedziale od 0 do 1, gdzie 0 oznacza brak odbicia (całe światło zostało pochłonięte), a 1 oznacza całkowite odbicie od powierzchni. W praktyce obiekty nie odbijają bądź pochłaniają stuprocentowo światła, niemniej sensory oraz procesy kalibracyjne nie są idealne, więc mogą pojawić się wartości odstające od tego przedziału, tak jak w tej sytuacji.

Można ten problem rozwiązać na dwa sposoby:

  1. Zastąpić te wartości brakiem danych (NA).
  2. Dociąć do minimalnej i maksymalnej wartości.

Pierwszy sposób może spowodować, że stracimy dużą część zbioru danych. Natomiast drugi sposób może powodować przekłamania.

# sposób 1
r = clamp(r, lower = 0, upper = 1, values = FALSE)
# sposób 2
r = clamp(r, lower = 0, upper = 1, values = TRUE)

Po przeskalowaniu wartości możemy wyświetlić kompozycję RGB. W tym przypadku zamiast funkcji plot() należy użyć funkcji plotRGB() oraz zdefiniować kolejność kanałów czerwonego, zielonego oraz niebieskiego. Często zdarza się, że kompozycje są zbyt ciemne/jasne, wtedy warto zastosować rozciągnięcie kolorów używając argumentu stretch = "lin" lub stretch = "hist".

# plotRGB(r, r = 3, g = 2, b = 1)
plotRGB(r, r = 3, g = 2, b = 1, stretch = "lin")

Klasteryzacja

library("cluster") # klasteryzacja danych

Dane do modelowania muszą zostać przygotowane w odpowiedni sposób. Modele klasyfikacyjne najczęściej na etapie trenowania wymagają macierzy lub ramki danych (data frame). Dane rastrowe można przetworzyć do macierzy przy użyciu funkcji values().

mat = values(r)
nrow(mat) # wyświetla liczbę wierszy
[1] 6000000

Za pomocą interaktywnej funkcji View() możemy sprawdzić jak wygląda nasza macierz.

View(mat)

W macierzy występują brakujące wartości. Zazwyczaj modele nie obsługują NA, więc musimy je usunąć. Służy do tego dedykowana funkcja na.omit().

mat_omit = na.omit(mat)
nrow(mat_omit)
[1] 5999675

Teraz przejdziemy do kolejnego etapu analizy, jakim jest wytrenowanie modelu. Istnieje wiele metod i modeli grupowania (patrz CRAN Task View), ale w tym przykładzie użyjemy prostego modelu grupowania metodą k-średnich. Ten model wymaga jedynie, aby podać z góry liczbę grup/klastrów (argument centers). Jest to algorytm stochastyczny, więc za każdym razem zwraca inne wyniki. Żeby analiza była powtarzalna musimy ustawić ziarno losowości – set.seed().

set.seed(123) # ziarno losowości
mdl = kmeans(mat_omit, centers = 3)

W wyniku powyższej operacji otrzymaliśmy m.in.:

  1. Obliczone średnie wartości grup dla poszczególnych kanałów (mdl$centers).
  2. Wektor ze sklasyfikowanymi wartościami macierzy (mdl$cluster).

Wyświetlmy te obiekty:

mdl$centers
        blue      green        red       nir
1 0.06072469 0.09172913 0.11470107 0.2175304
2 0.01790128 0.03401041 0.02827088 0.1998795
3 0.02791346 0.05896356 0.04148333 0.3397331
head(mdl$cluster)
[1] 2 2 2 2 2 2

Oznacza to, że pierwszy wiersz (reprezentujący pojedyncze oczko siatki) należy do grupy 2, drugi do grupy 2, trzeci do grupy 2, itd.

Walidacja

Nieodłącznym elementem modelowania jest walidacja opracowanych modeli. Wyzwaniem jest wybór właściwej metody grupowania dla konkretnego zbioru danych i określenie odpowiedniej liczby grup. Należy pamiętać, że zwiększenie liczby klastrów zwiększa podobieństwo między obiektami w klastrze, ale przy ich większej liczbie interpretacja staje się trudniejsza.

Najczęstszym sposobem walidacji wyników grupowania jest użycie wewnętrznych metryk, takich jak wskaźnik Dunna, wskaźnik Daviesa-Bouldina lub wskaźnik sylwetki (silhouette index). Na potrzebny niniejszej analizy użyjemy tego ostatniego.

Indeks sylwetki ocenia zbieżność i separację klastrów na podstawie odległości między obiektami w tym samym klastrze i między obiektami w różnych klastrach. Wartości tego wskaźnika mieszczą się w zakresie od -1 do 1. Wartość bliska 1 wskazuje, że obiekt jest dobrze zgrupowany i znajduje się daleko od sąsiednich klastrów. Wartość bliska -1 sugeruje, że obiekt mógł zostać przypisany do niewłaściwego klastra. Wartość bliska 0 wskazuje, że obiekt znajduje się bardzo blisko granicy pomiędzy różnymi klastrami. Ogólnie rzecz biorąc, wyższa wartość tego wskaźnika wskazuje na lepsze wyniki grupowania. Więcej szczegółów można znaleźć w dokumentacji cluster::silhouette().

Spróbujmy teraz obliczyć wartości tego wskaźnika. Zasadniczo wymaga to obliczenie podobieństwa każdego obiektu do każdego obiektu, co w naszym przypadku jest zadaniem niemożliwym (nasz zbiór danych składa się z ponad 6 mln obiektów). Aby to wykonać, musimy wykorzystać mniejszą próbkę (załóżmy \(n=10000\)). W funkcji musimy określić dwa obiekty –- wektor z klastrami oraz macierz niepodobieństwa, którą można wcześniej obliczyć za pomocą funkcji dist().

set.seed(123)
# losowanie indeksów
idx = sample(1:nrow(mat_omit), size = 10000)
head(idx)
[1] 2017743 3334670 4328362 3269750 2782437 5891675
# obliczenie wskaźnika sylwetki
sil = silhouette(mdl$cluster[idx], dist(mat_omit[idx, ]))
summary(sil)
Silhouette of 10000 units in 3 clusters from silhouette.default(x = mdl$cluster[idx], dist = dist(mat_omit[idx, ])) :
 Cluster sizes and average silhouette widths:
     3450      3423      3127 
0.4483002 0.4448400 0.4106674 
Individual silhouette widths:
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
-0.0601  0.3039  0.4940  0.4353  0.5926  0.6567 

Średnia zbieżność klastrów wynosi 0,44. Nie jest to najlepszy wynik (powinniśmy spróbować zwiększyć liczbę klastrów lub użyć innej metody grupowania). Możemy również zaprezentować wyniki na wykresie.

kolory = rainbow(3) # wybierz 3 kolory z wbudowanej palety `rainbow`
plot(sil, border = NA, col = kolory, main = "Silhouette Index")

Interpretacja

Istotą grupowania jest utworzenie grupy podobnych obiektów, natomiast naszym zadaniem jest zinterpretowanie tego, co reprezentują utworzone grupy i nadanie im nazwy. Interpretacja jest trudnym zadaniem i często wyniki są niejasne. W tym celu konieczne jest przeanalizowanie statystyk opisowych klastrów oraz wykorzystanie różnych wykresów i kompozycji map. Bardzo przydatna jest także znajomość właściwości spektralnych obiektów.

Spróbujmy zatem zinterpretować uzyskane skupienia, korzystając z wykresu pudełkowego. Największe możliwości wizualizacji danych dostarcza pakiet ggplot2. Tutaj można znaleźć darmowy podręcznik oraz gotowe “przepisy”. Wymieniony pakiet wymaga przygotowania zbioru danych do odpowiedniej postaci, tj. dane muszą być przedstawione jako ramka danych w tzw. formie długiej (wiele wierszy), podczas gdy standardowe funkcje do wizualizacji wymagają formy szerokiej (wiele kolumn). Takiej konwersji można dokonać w prosty sposób używając pakietu tidyr.

library("tidyr") # transformacja danych
library("ggplot2") # wizualizacja danych

Jak zauważyliśmy wcześniej, nasz zbiór danych jest dość duży i nie ma potrzeby prezentowania wszystkich danych. Możemy to zrobić efektywniej, używając wcześniej wylosowanej próby. Połączmy zatem wylosowane wiersze z macierzy z odpowiadającymi im klastrami (cbind()). Następnie macierz zamienimy na ramkę danych (as.data.frame()).

stats = cbind(mat_omit[idx, ], klaster = mdl$cluster[idx])
stats = as.data.frame(stats)
head(stats)
    blue  green    red    nir klaster
1 0.0728 0.0940 0.1100 0.2035       1
2 0.0121 0.0302 0.0217 0.3364       3
3 0.0258 0.0606 0.0341 0.3616       3
4 0.0166 0.0406 0.0253 0.3157       3
5 0.0770 0.1080 0.2156 0.2633       1
6 0.0368 0.0670 0.0572 0.3219       3

Wyświetlone dane mają formę szeroką (każdy kanał spektralny zapisany jest w osobnej kolumnie). Teraz musimy zmienić formę, w której otrzymamy dwie kolumn – kanał oraz wartość. W tym celu wykorzystamy funkcję pivot_longer().

stats = pivot_longer(stats, cols = 1:4, names_to = "kanal", values_to = "wartosc")

Dla formalności możemy jeszcze zmienić typ danych (klastrów i kanałów) na kategoryczny (factor). W praktyce związane jest to z uproszczeniem struktury danych (przejście ze skali ilorazowej do nominalnej).

stats$klaster = factor(stats$klaster)
stats$kanal = factor(stats$kanal)
head(stats)
# A tibble: 6 × 3
  klaster kanal wartosc
  <fct>   <fct>   <dbl>
1 1       blue   0.0728
2 1       green  0.094 
3 1       red    0.11  
4 1       nir    0.204 
5 3       blue   0.0121
6 3       green  0.0302

Ramka danych jest już przygotowana. Teraz stwórzmy prosty wykres pudełkowy.

ggplot(stats, aes(x = kanal, y = wartosc, fill = klaster)) +
  geom_boxplot()

Zmieńmy kilka domyślnych parametrów żeby poprawić odbiór ryciny.

etykiety = c("Niebieski", "Zielony", "Czerwony", "Bliska\npodczerwień")

ggplot(stats, aes(x = kanal, y = wartosc, fill = klaster)) +
  geom_boxplot(show.legend = FALSE) +
  scale_x_discrete(limits = kanaly, labels = etykiety) +
  scale_fill_manual(values = kolory) +
  facet_wrap(vars(klaster)) +
  xlab("Kanał") +
  ylab("Odbicie") +
  theme_light()

Na podstawie powyższego wykresu możemy przeanalizować właściwości spektralne klastrów, a tym samym zinterpretować, jakie obiekty reprezentują.

Finalna mapa

Ostatnim etapem jest stworzenie mapy klasyfikacyjnej pokrycia terenu na podstawie otrzymanego wektora z klastrami (mdl$cluster). Na początku musimy przygotować pusty wektor składający się z całkowitej liczby komórek rastra. Można to sprawdzić za pomocą funkcji ncell(). W naszym przypadku jest to 6 milionów komórek.

# przygotuj pusty wektor
wek = rep(NA, ncell(r))

Następnie musimy przypisać nasze grupy w wektorze w odpowiednie miejsca, tj. tym, które nie są zamaskowane (NA). Do niezamaskowanych wartości można odwołać się przez funkcję complete.cases().

# zastąp tylko te wartości, które nie są NA
wek[complete.cases(mat)] = mdl$cluster 

W ostatnim kroku należy skopiować metadane obiektu r, ale tylko z jedną warstwą, i przypisać mu wartości wektora wek.

# stwórz nowy raster
clustering = rast(r, nlyrs = 1, vals = wek)

Zaprezentujmy wynik grupowania na mapie, używając odpowiednich kolorów i nazw klastrów.

kolory = c("#d9d9d9", "#086209", "#2fbd2f")
kategorie = c("odkryta gleba", "lasy/woda", "roślinność")
plot(clustering, col = kolory, type = "classes", levels = kategorie,
     mar = c(3, 3, 3, 7))

Jeśli wynik jest zadowalający, możemy zapisać go za pomocą funkcji writeRaster(). Taki plik można później wczytać w R lub innym programie obsługującym dane przestrzenne (np. QGIS). Dodatkowo, w przypadku danych kategorycznych, podczas zapisu warto ustawić typ danych jako Byte / INT1U (pod warunkiem, że liczba kategorii nie przekracza 255).

writeRaster(clustering, "clustering.tif", datatype = "INT1U")

Zadanie

7. Pobierz scenę satelitarną o niskim zachmurzeniu i dotnij jej zasięg do wybranego powiatu. Wykonaj klasteryzację metodą kmeans oraz inną wybraną. Następnie dokonaj walidacji otrzymanych wyników wykorzystując wskaźnik silhouette. Przygotuj również wykres pudełkowy przedstawiający zmienność powstałych klastrów. Finalnie, zaprezentuj wynik klasteryzacji na mapie (dobierz odpowiedni schemat kolorów) oraz kompozycję RGB.